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Detailed protocol 


#Nichenet ligand-receptor interactions 
#This section covers Figure 6D and single cell RNA sequencing data deposited in GSE152042 


#Load packages 
library(Seurat) 
library(SeuratData) 
library(SeuratWrappers) 
library(cowplot) 
library(patchwork) 
library(ggplot2) 


library(nichenetr) 
library(tidyverse) 


#Prepare data - "epithelial" and "Seurat" objects correspond to the resclustering analyses performed in https://doi.org/10.7554/eLife.62810 
epithelial <- Renameldents(epithelial, '0' = "Basal 1", '1'= "Basal 2", '2' = "Inflammatory KC", '3' = "Inflammatory KC 2", '4' = "Differentiating KC", '5' = 
"Proliferating KC", '6' = "Inflammatory KC 3", '7' = "Inflammatory KC 4", '8' ="Proliferating KC 2", '9' ="Basal 3") 

Stromal <- Renameldents(Stromal, '0' ="Fb1", '1'="MyoFb", '2'="Fb2",'3'="Pericyte",'4'="Fb3", '5'="Fb4", '6' = "Fb5") 


DefaultAssay(Stromal) <- "RNA" 
DefaultAssay(epithelial) <- "RNA" 


#Rename conditions to "Healthy" and "Periodontitis" 
epithelial@meta.datal['new_ids']] <- apply(epithelial@meta.data, 1, function(x) {ifelse(x{['stim'] == 'Healthy’, 'Healthy’, 'Periodontitis')}) 
Stromal@meta.datal[[new_ids']] <- apply(Stromal@meta.data, 1, function(x) {ifelse(x[['stim'] == ‘Healthy’, 'Healthy’, 'Periodontitis')}) 


#Subset healthy stromal subtypes 

healthy_stromal <- subset(Stromal, idents = "Healthy’) 

DimPlot(healthy_stromal, reduction = "umap") 

Idents(healthy_stromal) <- "seurat_clusters" 

healthy_fibroblasts <- subset(healthy_stromal, idents = c("0", "2", "4", "5", "6")) 
Renameldents(healthy_fibroblasts, '0'="Fb1",'2'="Fb2",'4'="Fb3",'5'="Fb4",'6' = "Fb5") 


#Subset basal epithelial cells 
Basalt <- subset(epithelial, idents = "Basal 1") 


#The code below follows vignettes described in https://github.com/saeyslab/nichenetr: vignette("Seurat_wrapper", package="nichenetr"); 
vignette("seurat_steps", package="nichenetr") 


#Read in the NicheNet ligand-receptor prior network and ligand-target matrix 
ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds")) 
ligand_target_matrix[1:5,1:5] 


Ir_network = readRDS(url("https://zenodo.org/record/3260758/files/Ir_network.rds")) 
head(Ir_network) 


weiahted networks = readRDS(url("httos://zenodo.ora/record/3260758/files/weiahted networks.rds")) 


weighted_networks_Ir = weighted_networks$lr_sig %>% inner_join(Ir_network %>% distinct(from,to), by = c("from","to")) 
head(weighted_networks$lr_sig) 
head(weighted_networks$gr) 


#read expression data of interacting cells 
##seurat objects 
DefaultAssay(healthy_fibroblasts) <- "RNA" 
DefaultAssay(Basal1) <- "RNA" 


#Define a "sender/niche" cell population and a "receiver/target" cell population and determine which genes are expressed in both populations. In this study, the 
receiver cell population is the epithelial basal cells (Basal 1), whereas the sender populations are the stromal subpopulations. We hypothesised that Fb2 is a 
niche population. 


receiver = "Basal 1" 
expressed_genes_receiver = get_expressed_genes(receiver, epithelial, pct = 0.10) 
background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)] 


sender = "Fb2" 

expressed_genes_sender = get_expressed_genes(sender, healthy_fibroblasts, pct = 0.10) 
list_expressed_genes_sender = sender %>% unique() %>% lapply(get_expressed_genes, healthy_fibroblasts, 0.10) 
expressed_genes_sender = list_expressed_genes_sender %>% unlist() %>% unique() 


#Define a gene set of interest genes in the “receiver/target’ cell population that are potentially affected by ligands expressed by interacting cells (e.g. genes 
differentially expressed upon cell-cell interaction) 


condition_oi = "Periodontitis" 
condition_reference = "Healthy" 


DE_table_receiver = FindMarkers(object = seurat_obj_receiver, ident.1 = condition_oi, ident.2 = condition_reference, min.pct = 0.10) %>% 
rownames_to_column("gene") 


geneset_oi = DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_logFC) >= 0.25) %>% pull(gene) 
geneset_oi = geneset_oi %>% |[. %in% rownames(ligand_target_matrix)] 


#define a set of potential ligands which are expressed by the sender niche 
ligands = Ir_network %>% pull(from) %>% unique() 
receptors = Ir_network %>% pull(to) %>% unique() 


expressed_ligands = intersect(ligands,expressed_genes_ sender) 
expressed_receptors = intersect(receptors,expressed_genes_receiver) 


potential_ligands = Ir_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique() 


#perform NicheNet ligand activity analysis: rank the potential ligands 

ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = 
ligand_target_matrix, potential_ligands = potential_ligands) 

ligand_activities = ligand_activities %>% arrange(-pearson) %>% mutate(rank = rank(desc(pearson))) 

ligand_activities 


best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand) %>% unique() 


#check which sender population expresses these top-ranked ligands 
DotPlot(healthy_fibroblasts, features = best_upstream_ligands %>% rev(), cols = "RdYIBu") + RotatedAxis() 
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